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Abstract. Many stochastic differential equations that occur in financial modelling do 
not satisfy the standard assumptions made in convergence proofs of numerical schemes 
that are given in textbooks, i.e., their coefficients and the corresponding derivatives 
appearing in the proofs are not uniformly bounded and hence, in particular, not globally 
Lipschitz. Specific examples are the Heston and Cox-Ingersoll-Ross models with square 
root coefficients and the Ait-Sahalia model with rational coefficient functions. Simple 
examples show that, for example, the Euler-Maruyama scheme may not converge either 
in the strong or weak sense when the standard assumptions do not hold. Nevertheless, 
new convergence results have been obtained recently for many such models in financial 
mathematics. These are reviewed here. Although weak convergence is of traditional 
importance in financial mathematics with its emphasis on expectations of functionals 
of the solutions, strong convergence plays a crucial role in Multi Level Monte Carlo 
methods, so it and also pathwise convergence will be considered along with methods 
which preserve the positivity of the solutions. 



1. Introduction 
Consider the Ito stochastic differential equation (SDE) in M. d 

m 

dX t = a{X t )dt + bj{X t )dW t ij) , t e [0, T], X = x e R d (1) 

with drift and diffusion coefficients a, bj : M. d — > R d for j = l,...,m. Here Wt = 
(W^\ . . . , W^), t > 0, is an m- dimensional Brownian motion on a probability space 
(Q, J 7 , P) and superscripts in brackets label components of vectors. Throughout this 
article it will always be assumed that equation ([I]) has a unique strong solution. 

Explicit solutions of such equations are rarely known, thus one has to rely on numerical 
methods to simulate their sample paths X t (u>) or to estimate functionals E$(X) for some 
$ : C([0, T]; M. d ) — > R. Typically, such a numerical method relies on a discretization 

o < h < t 2 < . . . < t n = T 

and a global approximation on [0, T] is obtained by interpolation. 

In the case of the classical weak approximation the error of an approximation X to 
X is measured by the quantity 

|E0(X T ) - Ecj)(X T )\ 

for smooth functions : M d — > M. The test functions are a particular case of the 
general (path-dependent) functionals $. In the strong approximation problem the p-th 
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mean of the difference between X and X is analyzed, i.e. 

/ — \ i/p 

E sup \X tk -X tk \A 

v fc=0,...,n 7 

for the maximal error in the discretization points or 

/ , — , \ l f p 

E sup \X t -X t \ p ) 

v te[o,T\ J 

for the global error, where p > 1 and | • | denotes the Euclidean norm. Here the mean- 
square error, i.e. p = 2, is usually studied. The recent development of the Multi-level 
Monte Carlo method for SDEs [20] has revealed that strong error bounds are crucial 
for the efficient computation of functionals E$(X). 

While the strong error measures the error of the approximate sample paths X on 
average, the pathwise error is the random quantity 

sup \X tk (u) - X tk (u)\, u G Q 

k=0,...,n 

and 

sup \X t (cj) — X t (cu)\, uj G Q 
te[o,T] 

respectively. Here the error is analyzed for a fixed u G f2 without averaging. This 
quantity thus gives the error of the actually calculated approximation X tk (u), k = 
0, . . . , n, respectively X{u). 

The traditional weak and strong convergence analysis for numerical methods for sto- 
chastic differential equations (SDEs) relies on the global Lipschitz assumption, i.e. the 
SDE coefficients satisfy 

in 

\a{x)-a{y)\ + J2\ b j( x )- b M\ < L ■ \x - y\, x, y G R d 

for some L > 0. However, in many SDEs used for modelling in mathematical finance 
this assumption is violated, so the standard results (see [23111]) do not apply. 

The Constant Elasticity of Variance Model for asset prices [12] , which was introduced 
by Cox in 1975, is given by the SDE 

dS t = iiS t dt + a SI dW tl S = s > 

where fi G 1R, cr > and 7 G (0,1] and Wt, t > 0, is a one-dimensional Brownian motion. 
For 7 = 1 this is the standard Black-Scholes model (i.e. a geometric Brownian motion), 
while for 7 G (0, 1) the diffusion coefficient of this SDE is clearly not globally Lipschitz 
continuous. This SDE has a unique strong solution if and only if 7 G [1/2, 1] and takes 
values in [0, 00). 

The Ait-Sahalia model and its generalization [HHS], which are stochastic interest rate 
models, follow the dynamics 

dX t = (a^X^ 1 - a + a x X t - a 2 X r t )dt + aX p t dW u X = x > 

where ai,a,r,p > 0, i = — 1,...,2. Under certain conditions on the parameters (see 
[^5]). this SDE has a unique strong solution with values in (0, 00). Note that here the 
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diffusion coefficient grows superlinearly for large values of x while the drift coefficient 
has a singularity at x = 0. 

The Heston model [26] , which is an asset price model with stochastic volatility, is an- 
other example for an SDE with non-Lipschitz coefficients. This SDE takes non-negative 
values only and contains square root coefficients: 

dS t = fiS t dt + ^V t S t ( Vl-P 2 dW t {1) + p dW t {2) ) , S = s >0 

dV t = k(X - V t ) dt + 6y/V t dWl 2) , V = v > 0. 

The parameters satisfy p G R, k, A, 9 > and p G (—1, 1). The second component of 
this SDE is the Cox-Ingersoll-Ross process, which is also used as a short rate model |13j . 

Finally, the use of the inverse of the CIR process as volatility process leads to the 
so-called 3/2- model 

dS t = fiS t dt + ^/v t S t ( Vl-P 2 dW t {1) + p dW t {2) ) , S = s > 

dV t = Cl V t (c2 - V t ) dt + c 3 V? /2 dW t (2 \ Vo = v >0 

where ci, C2, c 3 > 0, see e.g. [27]. 

Motivated by these and other examples, the investigation of numerical methods for 
SDEs with non-Lipschitz coefficients has been an active field of research in recent years. 
This article, provides an overview of the new developments using the above equations 
as illustrative examples and discussing, in particular, Euler-type schemes. For some of 
the above equations exact simulation methods exist, see e.g. [HI [22] and also [9] for a 
class of one-dimensional equations, which are superior for the simulation of the SDEs 
at a single or a few time points. However, if a full sample path of the SDE has to be 
simulated or if the SDEs under consideration are part of a larger SDE system, then 
discretization schemes are typically more efficient. 



2. Pathwise Convergence Rates of the Euler Scheme and general 

Ito-Taylor Methods 

The pathwise error criteria are very robust with respect to the global Lipschitz as- 
sumption. One of the simplest approximation schemes for equation (JT]) is the Euler 
scheme 

m 

X tk+1 = X tk + a(X tk )A + b 3 (X tk )A k W«\ k = 0, 1, . . . , 

i=i 

with X = xo, where A = T/n, t k = kA and A k W = W tk+1 — W tk . The Euler scheme 
(and all other approximation methods that will be introduced below) depend on the 
stepsize A > 0, hence on n G N, but this dependence will be omitted whenever it is 
clear from the context. 

From the results of Gyongy [23] it follows that the Euler scheme has pathwise con- 
vergence order 1/2 — e also if the SDE coefficients are only locally Lipschitz continuous: 
for all s > 

sup \X tk -X tk \ <rjf -n^ 2+£ 

k=0,...,n 
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almost surely for a finite and non-negative random variable t]^ under the assumption 
that for all N G N there exist constants > such that 

m 

\a(x) - a(y)\ + ^ \ b j( x ) ~ b j(v)\ < L N ■ \x - y\, \x\, \y\ < N. 

3=1 

Thus, the pathwise convergence rate of the Euler scheme coincides up to an arbitrarily 
small e > with its strong convergence rate 1/2, but for the pathwise convergence rate 
no global Lipschitz assumption is required. 

Jentzen, Kloeden & Neuenkirch [37] observed that this is not a specific feature of the 
Euler scheme but, in fact, holds for general Ito- Taylor schemes of order 7 = 0.5,1.0,1.5,... 
For the definition of these schemes, see e.g. [39]. The Euler scheme corresponds to 
7 = 0.5, while 7 = 1.0 yields the Milstein scheme 

m m 

X tk+1 =X tk + a(X tk )A + J2bj(Xt k )A k WM+ £ LH h (X tk )I jlj2 (t k ,t k+l ) 

i=i 31,32=1 

with the differential operators 

L '=tt£i> j ' <<« 

fc=i 

and the iterated Ito-integrals 

I juj2 (s,t) = f f 2 dWW dW^\ j x ,j 2 = 1, . . . ,m. 

J s J s 

The Ito- Taylor scheme of order 1.5 is usually called the Wagner-Platen scheme. 

Theorem 2.1. Lei 7 = 0.5, 1.0, 1.5, .... Assume that a, b x , . . ., b m G C 2 ^ +1 (R d ; R d ) 
and moreover let X 1 ^ be the ltd- Taylor scheme of order'-/ with stepsize A = T/n. Then 
for every e > there exists a non-negative random variable rp e such that 



sup \X tk (u)-X t ' < r)](u) ■ n 

k=0,...,n 



IE 

-7+£ 



for almost all u G fi. 



The main ingredients to obtain this result are the Burkholder-Davis-Gundy inequality, 
which implies that all moments of an Ito-integral are equivalent, the following Borel- 
Cantelli-type Lemma, and a localization procedure. 

Lemma 2.2. (see [35],) Let a > 0, c p > for p > 1 and let (Z n ) neN be a sequence of 
random variables with 

(E\Z n \ p ) 1/P <c P -n- a 

for all p > 1 and n G N. Then for every e > there exists a finite and non-negative 
random variable rj e such that 

\Z n \ < f] e ■ n~ a+£ 

almost surely for all n G N. 
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The Burkholder-Davis-Gundy inequality and the Borel-Cantelli-type Lemma allow 
one to show that the Ito- Taylor scheme of order 7 has pathwise convergence rate 7 — e 
for smooth and bounded coefficients with bounded derivatives, thereby extending the 
classical mean-square convergence analysis in [39]. Then a localization argument is 
applied to avoid the boundedness assumptions. Roughly speaking, this localization 
argument works as follows: A fixed sample path X t (u), t G [0, T], of the SDE solution 
is bounded, i.e. stays in some open set B(u). However for the SDE 

m 

dY t = a(Y t ) dt + Y,h( Y t) dwP, Y = x 

with smooth and bounded coefficients a, bj with bounded derivatives, which coincide 
with the ones of the original SDE on B(u), the solution sample path Y t (cu),t G [0, T], 
coincides with X t (u),t G [0, T\. Asymptotically this also holds for the corresponding 
sample paths of the 7-Ito- Taylor schemes, so the pathwise convergence rates carry over. 

Note that all the examples of SDEs given in the introduction take non-negative val- 
ues only, so good approximation schemes should preserve this structural property. The 
(explicit) Euler scheme is, in general, not such a scheme, since its increments are condi- 
tionally Gaussian. For example, in case of the CIR process 

dX t = k(X -X t )dt + 9y/x' t dWt, X = x > 
the transition density of the Euler scheme reads as 

1 / (y - (x + k(\- x)A) 2 \ 

p(y; x) = = exp — - — , y G K, x > 0, 

V ; V2ne 2 xA \ 29 2 xA J' y 

so negative values can be obtained with positive probability even in the first step. This 
has lead to many ad-hoc corrections to prevent termination of the Euler scheme. The 
truncated Euler scheme 

x tk+1 = x tk + k(\ - x tk ) a + e^/xl A k W, k = 0, 1, . . . (2) 

was proposed in [2], while the scheme 

x tk+1 = x tk + k{\ - x tk )A + e^X^i A k W, k = 0, 1, . . . (3) 

was studied in [28J. Both approaches extend the mapping [0, 00) 3 x 1— > y/x G [0, 00) 
suitably to negative values of x. For the CIR process this idea was taken further by 
Lord, Koekkoek & van Dijk [30], who also proposed modifications of the drift coefficient 
for negative values of the state space. 

Example 2.3. The following table shows the average number of negative steps per 
path for the above Euler approximations of the CIR process. Scenario I (taken from 
[2]), corresponds to the parameters 

x = 0.05, k = 5.07, A = 0.0457, 9 = 0.48, T = 5 

while Scenario II (taken from (TTJ) uses 

x = 0.09, k = 2, A = 0.09, 9 = 1, T = 5. 



6 



PETER KLOEDEN AND ANDREAS NEUENKIRCH 



The stepsize for the Euler schemes is given by A = T/n with n = 512. 



average negative steps of / for 


Scenario I 


Scenario II 


Euler scheme (j2J) 


0.9141 


64.8611 


Euler scheme (j3J) 


1.0590 


74.5017 



The empirical frequency of negative paths is 0.4913 in Scenario I and 0.9990 in Scenario 
II. These results were obtained by a Monte Carlo simulation with iV = 10 6 repetition. 
They clearly indicate that the Euler scheme Q has a tendency for negative "excursions" . 
This can also be seen in Figured! which shows a sample path of the (linearly interpolated) 
Euler schemes ([2]) and ([3]) using the same path of the driving Brownian motion. The 
parameters used in this figure correspond to Scenario II. 

o 



For general SDEs the procedure of modifying the coefficients outside the support of 
the solution has been introduced systematically in [37]. For an SDE 

m 

dX t = a(X t ) dt + J2 h i ( x t) dW t (j) , X = x (4) 

3=1 

which takes values in a domain D C R d , i.e. 

P{X t e D, t > 0) = 1, (5) 

the auxiliary coefficients 

a(x) = a(x) ■ l D (x) + f(x) ■ l E (x), x E R d 

bj(x) = bj(x) ■ l D (x) + g (x) ■ l E (x), x G R d , j = l,...,m 

with E = R d \D are introduced there. A modified Ito-Taylor scheme of order 7 based 
on the auxiliary functions / and g is then the corresponding standard Ito-Taylor scheme 
for the SDE 



r{3) 

t > 



dX t = a{X t ) dt + J2 b j( X t) dWt ( 
3=1 

with a suitable definition of the derivatives of the coefficients on dD, see [37] for details. 
This method is well-defined as long as the coefficients of the equation are (27 + l)-times 
differentiable on D and the auxiliary functions are (27 — l)-times differentiable on E. 
The purpose of the auxiliary functions is twofold: to obtain a well-defined approximation 
scheme and to bring the numerical scheme back to D if it leaves D. In particular, the 
auxiliary functions can always be chosen to be affine or even constant. It was shown 
by Jentzen, Kloeden & Neuenkirch [37] that Theorem 12.11 adapts to modified Ito-Taylor 
schemes for SDEs on domains D <zM. d . 

Theorem 2.4. Let X be the solution of SDE (j4j) satisfying condition (J3]). Moreover let 
7 = 0.5, 1.0, 1.5, . . . and assume that 



a e C W {D- R d ), be C 2l+1 (D; m 
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Figure 1 . A path of Euler scheme fl2J vs. Euler scheme © for the CIR 
process and Scenario II 



and 



geC 2 ^\E- 



nd,m\ 



Finally let X 1,n be the modified Ito-Taylor method for X based on the auxiliary functions 
f and g with stepsize A = T/n. Then for every e > there exists a finite and non- 
negative random variable r\l: 9 e such that 



k=0,...,n 

for almost all u £ Q and all n£N. 



sup \X tk (oo)-X^(cu)\< V f : l 



u!) ■ n 



-7+£ 
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In the case of the Euler scheme, i.e. 7 = 0.5, the assumptions on a and b can be 
weakened to the assumption that a and b are locally Lipschitz continuous on D. For 
SDEs on domains in mathematical finance this condition is typically satisfied. In fact, 
in most cases the coefficients are infinitely differentiable. 

The CIR process satisfies 

P(X t > for all t > 0) = 1 

if and only if 2n\ > 9 2 . The latter assumption is typically satisfied in interest rate 
applications of the CIR process. Hence, modified Taylor schemes can be used here with 
D = (0,oo). The truncated Euler scheme ([2]) corresponds to the auxiliary functions 
f(x) = a(x), g(x) = 0, x < 0, while the scheme ([3]) uses the auxiliary functions f(x) = 
a(x), g(x) = a/— x, x < 0. Note that 2k\ > 6 2 is satisfied in Scenario I, but not in 
Scenario II. 

For structure preserving integration of the CIR process also the symmetrized Euler 
method 



k + l 



X tk + k(X - Xt k )A + 6JX tk A k W , k = 0, 1, . . . (6) 



was proposed in [TUl IB]. While the modified Euler schemes fl2]) and fl3]) may leave (0, 00) 
and are then forced back in the next steps, this scheme is always non-negative. Adapting 
this to general SDEs, which take values in a domain D, leads to the reflected Euler 
schemes, see e.g. jS], which are given by 

with Xq = xo, where 

m 

and a measurable projection function ifj : M. d \ D — > D U dD. A straightforward mod- 
ification of the above theorem yields a pathwise convergence order 1/2 — e for these 
reflected Euler schemes if the SDE coefficients are twice continuously differentiable on 
D. In the same way reflected Ito- Taylor schemes of arbitrary order can be constructed 
and analyzed. 

The symmetrized Euler scheme ([6]) corresponds to the reflection function ip(x) = \x\. 
The results on modified Ito- Taylor schemes and reflected Euler methods apply also to 
the generalized Ait-Sahalia model with D = (0, 00) if r > 1, p < (1 + r)/2, to the Heston 
model with D = (0, oo) 2 if 2k\ > 6 2 and to the 3/2-model with D = (0, oo) 2 and no 
further restrictions on the parameter. 

Example 2.5. To illustrate the above results consider Scenario I for the Cox-Ingersoll- 
Ross process. Figure [2] shows for two different sample paths u £ Q the maximum error 
in the discretization points, i.e. 

sup \X tk (u) -X tk (w)|, 

k=0,...,n 

Of 

(i) the truncated Euler scheme ()2]) 
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(ii) the symmetrized Euler scheme (EJ) 

(iii) the modified Milstein scheme with auxiliary functions f(x) = k(X 
i.e. a truncated Milstein scheme. 



■x), g(x) = 0, 



To estimate the pathwise maximum error for the above approximation schemes the Cox- 
Ingersoll-Ross process have been discretized with a very small step size using scheme (j2J). 
In Figure [2] log-log-coordinates are used, so the dots indicate the convergence orders 0.5 



10' 



10" 



10" 
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-trunc. Euler 
symm. Euler 
-trunc. Milstein 
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stepsize 
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FIGURE 2. Pathwise maximum error vs. step size for two sample paths 
for the Cox-Ingersoll-Ross model for Scenario I 



and 1. The pathwise convergence rates of all three approximation schemes are in good 
accordance with the theoretically predicted rates for moderate and small step sizes. For 
small step sizes both Euler schemes do not take negative values and hence coincide. 
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Moreover, for small step sizes the Milstein scheme is superior due to its first order 
convergence. 

o 

Numerical methods with pathwise convergence rates of high order are thus avail- 
able also for SDEs with non-globally Lipschitz coefficients. However, while pathwise 
convergence rates are very important for the analysis of random dynamical systems 
[3 ELY], one of the main objectives in mathematical finance is the pricing of (path- 
dependent) European-type derivatives, which means to compute real numbers E$(X) 
where $ : C([0,T];K ) — > R is the discounted payoff of the derivative. Since the in- 
tegrability of the random constants in the error bounds is an open problem, the above 
pathwise convergence rates do not imply weak or strong convergence rates. Neverthe- 
less, if <3> is bounded and continuous and if X 7 = (Xj) t e[o,T] is the piecewise linear 
interpolation of the 7-Ito- Taylor scheme (standard, modified or reflected) then 

E$(X 7 ) — ► E$(X) 

for n — > 00, so for bounded and continuous pay-offs (e.g. put options) one obtains 
at least the convergence of the corresponding standard Monte Carlo estimators for the 
option price. The same is true for barrier options with payoff of the form 

$(X) = (f}(X T )l{K 1 <\Xt\<K2, te[0,T]} 

with < Ki < K 2 < 00, if cj) is bounded and continuous and the law of sup tg j 0T ] \X t \ 
and inf ig [ 0j T] \X t \ has a density with respect to the Lebesgue measure. 



3. The Explicit Euler Scheme: Criteria for Weak and Strong 

Convergence 

It was shown by Higham, Mao & Stuart in [29] that the explicit Euler scheme 

fit 

X tk+1 = X tk + a(X tk )A + £ b 3 (X tk )A k W«\ k = 0, 1, ... , 

Xq = Xq 

is strongly convergent if the coefficients are locally Lipschitz continuous on M. d and a 
moment condition for the SDE and its Euler approximation is satisfied. This result can 
be extended to SDE on domains and the modified or reflected Euler scheme. 

Theorem 3.1. Let X be the solution of SDE (jl]) satisfying condition ([5]). Moreover, 
let X n be the modified Euler scheme based on the auxiliary functions / G C , (£';lR d ), 
g G C(E; M. d ' m ) with stepsize A = T/n or let X n be the reflected Euler scheme based on 
the projection function ^:£->DU dD with stepsize A = T/n. Assume that 

aeC 2 {D-m d ), b eC 2 (D-m d > m ) 

and furthermore, assume that for some p > 2 

supE max \X^ \ P + E sup \X t \ p < 00. (7) 

ngN k=0,...n te[0,T] 
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Then 

lim E 

k=0,...,n 



lim E max \X tk - X™ 



Proof. From the results of the previous section 

lim max \X f , — X? I = 

n— >oo k=0,...,n 

hold, almost surely. However, assumption (J7|) implies the uniform integrability of 

max \X t t -X n | 2 , n G N. 

The assertion now follows, since uniform integrability allows integration to the limit. □ 

Note that assumption (JTj) is easily verified if the SDE coefficients have linear growth 
on D, i.e. 

m 

H*)\ + ^2\bj(x)\ < C- (1 + |x|), xeD, 

3=1 

for some C > 0. Turning back to the Cox- Ingersoll- Ross process this gives us strong 
convergence of the Euler schemes (j2J), ([3]) and under the assumption 2k\ > 2 . Note 
that for the Euler schemes ([2]) and strong convergence without a restriction on the 
parameter has been shown in [2] and [2B] using a Yamada function technique. This 
technique has also been applied by Gyongy & Rasonyi in [21] to obtain the following 
result: 

Theorem 3.2. Let a±, a 2 , b : R — > R. Consider the one-dimensional SDE 

dX t = {axiXt) + a 2 (X t )) dt + b(X t ) dW u t G [0, T], I = i 6l 

and let X™ be the corresponding Euler scheme with stepsize A = T/n. Moreover, let a 2 be 
monotonically decreasing and assume that there exists constants a G [0, 1/2], (3 G (0, 1] 
and C > such that 

|ai(x) - ai(j/)| < C • |s - |a 2 (x) - a 2 (j/)| < C- |s - yf, 

\b(x)-b{y)\ < C- \x-y\^ +a 
for all x, y G R. Then, for all p G N, i/iere exzsi constants K^ 13 > sitc/j that 

Ky ■ — ^ /or a = 



E , m n aX „ ^ - P ^ 1 ^ • + M f° T a G (°> V2) 

^ " + M « = 1/2 



This result can be applied to the CEV model 

dX t = fiX t dt + aX? dW u 

if the mapping [0, oo) 9 i 4 i 7 6 [0, oo) is extended to (— oo,0), e.g. as (x + ) 7 or \x\ 
Theorem 13.21 then yields strong convergence of the corresponding Euler schemes. 
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Example 3.3. Whether the convergence rates predicted from Theorem 13.21 are sharp 
for the CEV model remains an open problem. The following simulation study suggests 
that the Euler scheme has strong convergence order 1/2, at least for some parameter 
constellations. To better preserve the positivity of the CEV process, the Euler scheme 
is applied to the SDE 

dX t = fj,\X t \ dt + a{X+y dW t 

which still fulfills the assumptions of Theorem 13.21 with ai = 0, i.e. (3 = 1, and a = 
7 — 1/2. Its solution coincides with the CEV process. 

Figure [3] shows the empirical root mean square maximum error in the discretization 
points versus the step size for the parameters 

Set I: /1 = 0.1, a = 0.3, 7 = 0.75, T = 1, x = 0.2 
Set II: fi = 0.2, a = 0.5, 7 = 0.55, T = 1, x = 0.5 

The empirical mean square maximum error in the discretization points is estimated by 

/ TV \ V2 

— > max A. — a, 

I N^k=0,...,n l tk k J 

with N = 5 • 10 4 . Here X* is the numerical reference solution obtained by using the 

same Euler scheme with very small step size and X*'( t \~X n ' ^ are independent copies of 
X*,X . For both sets of parameter a good accordance with the convergence order 1/2 
is obtained. (The dots in the figure indicate convergence order 1/2). A regression of 



Euler scheme - set I 
-Euler scheme -set II 



1(f - 



10- 3 10- 2 



stepsize 

Figure 3. Root mean square error of the Euler scheme vs. step size for 
the CEV process for the parameter sets I and II 

the numerical data yields moreover the empirical convergence order 0.493923 for set I, 
respectively 0.509903 for set II. 

o 
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But do Theorems 13.11 and 13.21 have any consequences for the other examples? Unfor- 
tunately not: for the Heston, Ait-Sahalia and 3/2-models, no linear growth condition is 
satisfied. Even worse, for the Ait-Sahalia model and the 3/2-model the moments of the 
Euler scheme explode! In the case of the 3/2-model the latter can be deduced from the 
following Theorem, which was obtained by Hutzenthaler, Jentzen & Kloeden in [31]. 

Theorem 3.4. Let a, b : R — > R and assume that the one- dimensional SDE 

dX t = a(X t ) dt + b(X t ) dW u t e [0, T], X = x e R 
has a unique strong solution with 

sup E\X t \ p < oo 
te[o,T] 

for one p G [1, oo). Moreover, let b(xo) ^ and let C > 1, /3 > a > 1 be constants such 
that 

max( \a(x)\ , \b(x)\ ) > — • \xf and min( \a(x)\ , ) < C ■ \x\ a 

for all \x\ > C. Then, the corresponding Euler scheme X" with stepsize A = T/n 
satisfies 

lim E|X r - T T \ V = oo and lim |E|A: T | P - E|A™ | p | = oo. (8) 
In the case of the 3/2-model, which has finite moments up to order p < 2 + the 

c 3 

coefficients are 

a(x) = —c 2 x 2 + Cic 2 x, b(x) = c 3 (x + ) 3 ^ 2 , x G R, 

so the assumptions of the above Theorem are satisfied for a = 3/2, (3 = 2 and C 
sufficiently large. 

Concerning the Ait-Sahalia model, the moments of the Euler scheme already explode 
in the second step. Here the first step of the Euler scheme has a Gaussian distribution 
with mean x + («-i^o 1 — «o + a\XQ — ot2X r ) A and variance a^x^fA. The inverse of the 
first step must be computed for the second step of the Euler scheme, so the moments 
of the second step are infinite, since inverse moments of a Gaussian random variable do 
not exist. 

Why the moments of the Euler scheme diverge for superlinearly growing coefficients 
- even without a singularity - can be nicely illustrated by considering the SDE 

dX t = -X? dt + adW t , X = x (9) 

with a > for which the Euler scheme reads as 

X k+1 =X k {^-\X k \ 2 ^)+^ k w. (io) 

In the deterministic case, i.e., fl9]) and ( TlT)]) with o = 0, the Euler approximation of 
the deterministic equation is known to be unstable if the initial value is large (see e.g. 
Chapter 6 in [IS]). For example, if Xq = n, T = 1 and A = n~ l then 

*" h = »(i-£)—' 

and therefore 

x£=^(l-K| 2 A)«n 8 . 
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Iterating this further, one obtains 

>J 2fc ) 

for k = 0, 1, . . . ,n. Thus, X™ n grows double-exponentially fast in n. In the presence of 
noise (a > 0) there is an exponentially small event that the Brownian motion leaves the 
interval [— 2n, 2n] and on this event the approximations grow double-exponentially fast 
due to the deterministic dynamics. Consequently this double-exponentially growth can 
not be compensated by the exponentially small probability of this event, which leads to 
the moment explosion of the Euler approximation. 

Example 3.5. That rare events lead to the explosion of the moments of the Euler 
scheme can be also seen from the following numerical example. Consider the volatility 
process in the 3/2-model 

dV t = Cl V t (c 2 - V t ) dt + c 3 V t 3/2 dW t , V = v >0 

with 

c x = 1.2, c 2 = 0.8, c 3 = 1, T = 4, v = 0.5 
and try to compute 

E\X T \ = 0.566217 
using the standard Monte Carlo estimator 

1 N 

i=l 

— n.(l) — n,(N) — n 

where X T ' , . . . ,X T ' are iid copies of X T . The exact value for E|X^| is computed 
using the inverse moments of the CIR process, see e.g. [21]. While for a moderate 
number of repetitions the estimator seems to converge for small step sizes (and the 
'Inf '-outputs seem to be some numerical instabilities due to the large step sizes), the 
estimator explodes even for small step sizes when increasing the number of repetitions 
- as predicted by Theorem 13.41 Despite of this the Euler scheme for this SDE converges 
pathwise with rate 1/2 — e due to Theorem 12.41 



Repetitions iV / stepsize A 


2° 


2" 2 


2 -4 


2-6 


2- 8 


2 -io 


10 3 


6.327232 


Inf 


Inf 


0.550185 


0.553499 


0.555069 


10 4 


6.894698 


Inf 


Inf 


Inf 


0.562716 


0.563352 


10 5 


7.430606 


Inf 


Inf 


Inf 


0.566218 


0.567106 


10 e 


7.227379 


Inf 


Inf 


Inf 


Inf 


0.565750 


10 Y 


7.279187 


Inf 


Inf 


Inf 


Inf 


Inf 



A similar moment explosion arises if a Multi-level Monte Carlo method is used to 
estimate E|Xy|. This is shown, also for more general SDEs, in [3"E] . 

o 



However, in some cases using the Euler scheme one still obtains a convergent Monte 
Carlo estimator for functionals of the type E0(Xr). The standard Euler-based estimator 
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for the latter quantity is 



N 



In the classical case, i.e. if a, b, <ft e C 4 (1R; 
tives and a, b globally Lipschitz, one has 



i N 



(ii) 



with at most polynomially growing deriva- 



E 



1 

N 



N 



(T/) - E0(X 



< if 



Bias 



+ K 



MC 



1 

iV' 



see e.g. |39j . The first term on the right hand side corresponds to the squared bias of 
the Euler scheme, while the second term corresponds to the variance of the Monte Carlo 



simulation. It is thus optimal to choose iV 



ir 



for balancing both terms with respect 



to the computational cost (number of arithmetic operations, function evaluations and 
random numbers used), see [IS]. The corresponding Monte Carlo estimator has then 
convergence order 1/3 in terms of the computational cost. 

Hutzenthaler & Jentzen could show in [32J that if the global Lipschitz assumption on 
the drift-coefficient is weakened to 



- y){a(x) - a(y)) <L(x — y)' 



x,y e 



(12) 



for some L > 0, then one still has 



N 



almost surely for all e > and almost-surely finite and non-negative random variables 

Ve- 

Weak approximation under non-standard assumptions is also studied by Milstein & 
Tretyakov in [32]. In their approach, simulations which leave a ball with sufficiently 
large radius are discarded. In the context of the Euler scheme with equidistant stepsize 
this estimator reads as 



1 N 



|X"' W |<R}- 



For coefficients a, b and functions <fi satisfying a Lyapunov-type condition still a con- 
vergent Monte Carlo estimator is obtained, when matching the discarding radius R 
appropriately to the number of repetitions N and the stepsize of the discretization n. 

Condition (jl2j) on the drift coefficient is the so-called one-sided Lipschitz condition. 
This condition is also very useful to obtain strong convergence results for implicit Euler 
methods and tamed Euler schemes, which will be explained in the next section. Very 
recently a unifying framework for the analysis of Euler-type methods has been provided 
in 
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4. Strong convergence of implicit and tamed Euler schemes 



The condition in Theorem 13.11 for the strong convergence of the Euler scheme which 
is usually difficult to verify is the finiteness of its moments, i.e. 

sup E max \X" ' \ p < oo 

„ £ N k=0,...,n 

for some p > 2. Moreover, this condition may even fail to hold for specific equations, see 
Theorem 13 .41 However, both problems can be overcome in some situations if appropriate 
drift-implicit Euler schemes are used. The split-step backward Euler scheme is defined 
as 

m 

x* = x tk + a(x; k )A, X tk+1 = X* + £ b^X^A.W^ (13) 

for k = 0, 1, . . . with X = xq, while the backward or drift-implicit Euler scheme reads 
as 

m 

X tk+1 = X tk + a(X tk+1 )A + £ bj(X tk )A k WW. (14) 

j'=i 

Both schemes are defined via an implicit equation, whose solvability relies on the prop- 
erties of the drift-coefficient a. The following result has been obtained by Higham, Mao 
& Stuart in [29] . 

Theorem 4.1. Let a,bj G C^R"*; R d ), j = 1, . . . ,m, and assume that there exist con- 
stants Li, L 2 > such that 

(x - y, a(x) - a{y)) < L x ■ \x - y\ 2 , x,y E R d , 

m 

\bj(x) - b 3 (y)\ 2 <L 2 -\x- y\ 2 , x,y e R d . 

i=i 

Then, the split-step backward Euler scheme given by ( fT3|) with stepsize A = T/n is well 
defined for A < A* := 1/ max{l + 2L±, AL2} and satisfies 

lim E max \X tk - XT. I 2 = 0. 



n— >oo k=0,...,n 



The conditions on the coefficients imply that the SDE has bounded moments of any 
order, and also allow one to show that the split-step Euler method has moments of 
any order. The implicitness of the method is crucial for the latter. Furthermore, the 
split-step Euler method coincides with the explicit Euler method for the perturbed SDE 

m 

dX t A = a{h A (X t A )) dt + J2 bj{h A (X A )) dW {j \t), X A = x . (15) 

Here the function h\ : R d -> R d is defined as the unique solution of the equation 

h A (x) = x + a(h A (x))A, x G R d 

with A < A*. Since h A converges to the identity for A — > this perturbed SDE is close 
to original SDE. To establish Theorem I4.1[ it thus remains to show that the split-step 
backward Euler scheme is close to (fl~5|) . which can be done along the lines of the proof 
of Theorem 13.11 
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If the drift-coefficient is additionally also polynomially Lipschitz, then the standard 
strong convergence rate 1/2 can even be recovered. 



Theorem 4.2. Let the assumptions of Theorem 4-1 hold and assume additionally that 
there exist C, q > such that 

\a(x) - a{y)\ < C • (1 + \x\ q + \y\ g ) -\x-y\, x,y e M. d . 

Then, the split-step backward Euler scheme given by ffT3l) and the backward Euler scheme 
given by ()14|) are well defined for A < A* and have strong convergence order 1/2, i.e. 
for both schemes there exists a constant K > such that 

X,, - x n 

k=0,...,n 



E max \X U - XI I 2 < K ■ n" 1 . 



As pointed out above, in each step of both schemes an implicit equation has to be 
solved. If the function h& is not known explicitly, this has to be done numerically and 
may be time-consuming. Solving implicit equations can be avoided by using the so-called 
tamed Euler method, which has been proposed by Hutzenthaler, Jentzen & Kloeden in 

^ = ^ + ixi fir MA a ^ A + E bj(X tk )A k W^. (16) 
l + \a{A tk )\a j=1 

Here the drift-term is "tamed" by the factor - — .4 . . . in the k-th step, which prevents 

J l+|o(X tj .)|A ^' 

a possible explosion of the scheme. 

Theorem 4.3. Let the assumptions of Theorem \4-S\ hold. Then, there exists a constant 
K > such that the tamed Euler scheme given by (fl6l) satisfies 

\X f , -x n 

k=0,...,n 



E max \X th - X? I 2 < K ■ n" 1 . 



Here, the difficulty is again to control the moments of the approximation scheme. For 
this appropriate processes are used that dominate the tamed Euler scheme on subevents 
whose probabilities converge sufficiently fast to one. 

The Theorems given so far in this section require the diffusion coefficient to be glob- 
ally Lipschitz, which is often not fulfilled in SDEs arising from mathematical finance. 
However, the backward Euler method can be also successfully applied to the Ait-Sahalia 
interest rate model 

dX t = (a^X" 1 -a + a x X t - a 2 X r t ) dt + aX p t dW t (17) 

where at, a > 0, % = — 1, . . . , 2 and r, p > 1. The following result has been obtained by 
Szpruch et al. in 



Theorem 4.4. Consider the SDE (1171) and assume that 

r + 1 > 2p. 

Then the corresponding backward Euler method ( 1141) with stepsize A = T/n is well 
defined if A < \jct\, and 



lim E max \X U - X" I 2 = 0. 

k=0,...,n " k 
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Here the drift coefficient is still one-sided Lipschitz on the domain of the SDE, i.e. 
(x - y)(a(x) - a(y)) < a x \x - y\ 2 , x, y > 0, 
and, moreover, —a is coercive on (0, oo), i.e. 

lim a(x) = oo lim a(x) = — oo. 

x— x—too 

These two properties ensure that the drift-implicit Euler scheme for 017p is well-defined 
and, in particular, takes only strictly positive values. 

The drift coefficient in the volatility process 

dV t = Cl V t (c 2 - V t ) dt + c 3 V t 3/2 dW tl V = v >0 

in the 3/2-model is also one-sided Lipschitz on (0, oo). It does not, however, satisfy the 
coercivity assumption. Consequently, the drift-implicit Euler scheme cannot be applied 
here, since the implicit equation may not be solvable. Note that very recently, Higham 
et al. introduced in [30] a double-implicit Milstein scheme, which is strongly convergent 
for the 3/2-model and similar SDEs. 



5. Strong Convergence Rates for the approximation of the 
Cox-Ingersoll-Ross process and the Heston model 

Strong convergence rates for the approximation of the CIR process 

dX t = k(\ - X t ) dt + 9^/X~ t dW t , t G [0, T], X = x > 

with k, A, 9 > have been a long standing open problem, even in the regime where the 
CIR process does not hit zero, i.e. when 2k\ > 9 2 . 

The first non-logarithmic rates were derived by Berkaoui, Bossy & Diop for the sym- 
metrized Euler scheme (El), i.e. 



X tk + K(X-X tk )A + 9 x X tk A k W 



They showed in [8] that 



E max \X ft - X fl \ 2p < C„ ■ A p 

k=0,...,n 



under the assumption 

2k\ r- [ \[k 



> 1 + V8 max WB-y/l&p-l, 16p - 2 j 



ff 2 

where the constant C p > depends only on p, k, A, 9, Xq and T. Strong convergence rates 
for a drift-implicit Euler-type scheme were recently obtained under mild assumptions by 
Dereich, Neuenkirch & Szpruch in [15]. Their key tool is the use of the Lamperti- 
transformation: by the Ito formula, the transformed process Y t = yXt satisfies the 
SDE 

a 

dY t = —dt + (3Y t dt + ~fdW t , t > 0, Y = ^ (18) 
Yt 

with 

4k\-9 2 k 9 

a = , p = , 7 = -. 

8 2 ' 2 
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At first glance this transformation does not help at all, since the drift coefficient of the 
arising SDE is singular. However, 

cx 

a(x) = — h fix, x > 0, 
x 

satisfies for a > 0, G R the restricted one-sided Lipschitz condition 

(x - y)(a(x) - a(y)) < /3(x - y) 2 , x,y > 
The drift-implicit Euler method with stepsize A > in this case is 

Y tk+1 =Y tk +(=?- + (3Y tk+1 ) A + 1 A k W, k = 0, 1, 
fc+i 



with Y = Jx~$, which has the explicit solution 



Y tk + 1 A k W , (Y tk + 1 A k W) 2 , «A 



2(1 -/3A) V 4(1 -/?A) 2 1-/3A' 



Setting 



X, fc =F, fc , A; = 0,1,..., (19) 

gives a positivity preserving approximation of the CIR process, which is called drift- 
implicit square-root Euler method. This scheme had already been proposed in [1], but 
without a convergence analysis. Piecewise linear interpolation, i.e. 

X t = — X tk H -^X tk+1 , t G 

gives a global approximation (X t )tG[o,T] of the CIR process on [0, T]. The main result of 
is: 



Theorem 5.1. Let 2n\ > 9 2 , x > and T > 0. Then, for all 

1 < p < 



there exists a constant K p > such that 



2 

i/p 



(e max \X t - X t A ? < K p ■ ^/fog(Kj\ ■ VA, 
/or all A E (0,1/2]. 

The restriction on p arises in the proof of the convergence rate when controlling the 
inverse p-th moments of the CIR process, which are infinite for p > 2kX/8 2 . For further 
details, see [15]. Note that for SDEs with Lipschitz coefficients the convergence rate 
\J\ log(A)| • \/A is best possible with respect to the above global error criterion, see 



So the convergence rate given in Theorem 15.11 matches the rate that is optimal under 
standard assumptions. 

Other approximation schemes for the strong approximation of the CIR process can 
be found in [U |2U [2S] • Among them is the drift-implicit Milstein scheme 

Z tk+1 = Z tk + «(A - Z tk+1 )A + 9\J~2t k A k W + j ((A k W) 2 - A) 
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with Zq = Xq, see [2T]. It can be rewritten as 

= TT^a (v^ + W + TT^a ( KA - 7 ) A ' < 20) 



so this scheme preserves the positivity of the CIR process if AkX > 6 2 . It coincides up to 
a term of second order with the drift-implicit square-root Euler method, since the latter 
can be written as 

= l —r ( + -A k w] + ^— (kX - — ) A 

fe+1 1 + kA \ V tk 2 k J 1 + kA\ 4 

1 + kA \ o 2 v tfc+1 



Moreover the drift-implicit Milstein scheme dominates the drift-implicit square-root Eu- 
ler method: 

Lemma 5.2. Let 2k\ > 9 2 , x > and T > 0. Then 

P(Z tk >X tk , A; = 0,1,...) = 1. 
Proof. The numerical flow for the drift-implicit Milstein scheme is given by 

^ k ' A > = 17727 + l A ' w ) 2 + 177a ( kA " t) A 

and for the drift-implicit square-root Euler method it satisfies 



2 

Wy(i; k, A) H = v <p-v(a;; k, A) ) A 

1 + ^A I 8y/tpx(x;k,A) 2 V 



2 



^ + ~A k w)\ — i-j- f«A - ^ ) A. 



1 + ftA \ v 2 J 1 + kAV 4 

From [1] it is known that ipx is monotone, i.e. 

(Px(xx; k, A) > (fx(x 2 ; k, A) 

for xi > X2- Thus it remains to show that 

(Pz(x; k, A) > ip^x; k, A) 

for arbitrary A > 0, k — 0, 1, . . ., x > 0. However, this follows directly by comparing 
both flows. □ 

The above property allows one to show the strong convergence of the drift-implicit 
Milstein scheme, which seems not to have been established yet in the literature. 

Proposition 5.3. Let 2k\ > 9 2 , x > and T > 0. Then 

lim E max \X t . -Z™\ 2 = 0. 

n— too k=0,...,n k 
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Proof. First note that the drift-implicit square-root Euler method can be rearranged as 



X tk+1 = Vz{X tk - k, A) - K [ h+ \x t - X tk+1 ) dt 







■k+1 



ff 2 



(y/Xt - y/X th )dW t - -(A k W 2 - A) 



where is the numerical flow of the drift-implicit Milstein scheme defined in the proof 
of the above Lemma. Thus the error e k = X tk — Z tk satisfies the recursion 



e k+1 = e k - Ke k+i A + d(y/X^ k - \J~Z^\A k W + p k+1 
with eo = 0, where 



(21) 



p fc+ l = -K 

Now (J2U) gives 



tk+l 



(X s - X tk+1 ) ds + 9 



tk+i 



'X s -y/X tk )dW s . 



1 + kA 



e k + 6[ JX tk - \ Z tk )A k W + p k+1 



so 



pa- 



e=o K ' i=o v 



A) 



Straightforward calculations using ( 120]) yield 

sup sup EZ tfc < oo. 

nSN fc=0,...,n 

Then applying the Burkholder-Davis-Gundy inequality to the martingale 

fc-i 

M k = ^(1 + KA) e (y/xT t ~ yjZt^AtW, k = 0, 1, . . . 



1=0 



gives 



n-1 



E sup e 2 k < c^(1 + kA) 2£ E 



k=0,...,n 



e=o 



Xt e — y Z te 



A 



(22) 



+ c E sup 

k=l,...,n 



k-1 



5 a + «A)"-' Pm 



£=0 



k-V 



Here and below constants whose particular value is not important will be denoted by c 
regardless of their value. 

It remains to estimate the terms on the right side of the equation f[2"2"j) . The previous 
Lemma implies that 



E l^t fc _ x tk \ - E(Z tfc - x 



tk/1 



so 



nZ tk - X tk \ < 2 nx tk - X tk \ + \E{Z tk - X tk ] 
Clearly, Theorem 15.11 yields 

max E\X tk -X tk \<c - y/\ log(A)| • VA. 

k=0,...,n 
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Moreover, 

EZ tk+1 =EZ tk + K(X-EZ tk+1 )A 1 

which is the drift-implicit Euler approximation of 

-t 



EX t = x + [ k(X- EX s ) ds, t e [0, T], 
Jo 



at the discretization points = kA, so 



max \E{X tk -Z tk )\<c-A. 

fc=0,...,n 



Hence 



which gives 



mx -E\X tk - X tk | < c • VIMA)! • Va 



+ KA) 2k E 



fc=0 



A<c- v /|log(A)|-V / A 



(23) 



since \^fx — ^/y\ < yj\x — y\ for x, y > and sup neN sup fc=0 n (l + nA) 2k < oo. 

For the second term, applying the Burkholder-Davies-Gundy inequality and Jensen's 
inequality yield 



E sup 

k=l....,n 



fc-1 



S (1 + kA)^ P£+1 



n-l 



fc=0 



Now 

so it follows that 



tk 



(X t -X tk+1 )dt 
E\X t -X s \ 2 <c-\t-s 



2 n-l 
+ C 

fc=0 



tk + l 



E E 

=0 

s,te [o,T] 



dt. 



E sup 

fc=l,...,n 



fc-1 

^ (1 + K, 



A) 



k-i 



Pt+1 



Cc-VA, 



which completes the proof of the proposition. 



(24) 
□ 



Alternatively, Proposition 15.31 could have been obtained by deriving the pathwise con- 
vergence of the drift-implicit Milstein scheme and establishing the uniform integrability 
of the squared maximum error. Note that the above proof gives also the convergence 
order 1/4 up to a logarithmic term. However this rate seems to be suboptimal, see the 
following numerical example. 

Example 5.4. The Figures H] and |5] show the empirical root mean square maximum 
error in the discretization points, i.e. 

1/2 



if; 

N ^k 

i=l 



fc=0,...,n 



tfc ife 



versus the step size for the approximation of the CIR process. Consider the 
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(i) truncated Euler scheme (j2J) 

(ii) drift-implicit square-root Euler (fl9|) 

(iii) drift-implicit Milstein scheme (120]) 

for the Scenarios I and II (see Example [23]). Scenario I satisfies the condition of Theorem 
15. H and Proposition l5.3l since 2k\/8 2 = 2.011276 . . ..This condition is violated in Scenario 
II where 2k\/6 2 = 0.36. In the latter scenario, the truncation \/x+ in the definition of 
the schemes (ii) and (iii) is used, since both discretization schemes may take negative 
values here. The numerical reference solution X* is computed in Scenario I using scheme 
(ii) with very small stepsize and in Scenario II with scheme (i) with a very small stepsize. 
The number of repetitions of the Monte Carlo simulation is N = 5 • 10 4 . In the log- 




FlGURE 4. Root mean square errors vs. step size for Scenario I of the CIR process 

log coordinates here, the dots indicate the convergence orders 0.5 and 1.0 in Figure H] 
and 0.25 and 0.5 in Figure [5] respectively. For Scenario I the empirical mean square 
error for the truncated Euler scheme seems to decay with the order 0.5, while the other 
schemes seem to have an empirical convergence order close to 1.0. (For smooth and 
Lipschitz coefficients the Milstein scheme is of order one for the maximum error in the 
discretization points.) 

For Scenario II, these convergence orders deteriorate and for all schemes are signifi- 
cantly lower than one half, see also the following table, where the convergence orders 
have been estimated by a linear regression. 



empirical conv. order / for 


Sc. I 'part' 


Sc. II 'part' 


Sc. I 'full' 


Sc. II 'full' 


truncated Euler 


0.5739 


0.3193 


0.6446 


0.2960 


drift-imp. square-root Euler 


0.9281 


0.2734 


0.8491 


0.2837 


drift-imp. Milstein 


0.9447 


0.3096 


0.8719 


0.2871 



Here 'part' denotes the results for the linear regression using only the step sizes A = 
5 • 2~ J , j = 7, ... ,13, while 'full' uses the full data set, i.e. the step sizes A = 5 ■ 2~- J , 
J = 4, . . . , 13. 
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— trunc. Euler 

— imp. Milstein 
— sqr. Euler 



IIP 
stepsize 



Figure 5. Root mean square errors vs. step size for Scenario II of the CIR process 



Applying the Lamperti-transformation also to the asset price in the Heston model 
gives the log-Heston model 



d\og(S t ) =L- ±Y?j dt + Y t (Vl - P 2 dW t {1) + pdW^) 



dY 



AkX-6 2 1 k\ , 9 7TjA2) 



t 1 



S = s > 
Y = y > 0. 



The approximation of the log-Heston price is then a simple integration problem. Using 
the Euler scheme for the log-price equation and the drift-implicit square-root Euler 
scheme for the volatility process yields an approximation H tk of log (St fc ) given by 

k— 1 / -i \ k—1 
£=0 ^ ' 1=0 

This is extended by piecewise linear interpolation to [0, T]. 

Corollary 5.5. Let 2k\ > 9 2 , x > and T > 0. Then, for all 

2k\ 



there exists a constant K p > such that 



(e max | \og(S t ) - H t A * < K p ■ V\hg(A)\ ■ VX, 
V te[o,T] J 

for all A E (0,1/2]. 
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Note that in the Heston model moment explosions may appear according to the pa- 
rameters of the SDE. In particular, for p > 1 one has ESf < oo for allt > if and only 
if 

For more details see e.g. [6]. Whether this phenomenon also arises for discretization 
schemes for the Heston model is unknown at the time of writing. 

Example 5.6. In this example we test the efficiency of the Multi-level Monte Carlo 
estimator P m i see [TH1 ED], based on the above approximation scheme for the valuation 
of a European Call option, i.e. for 

p = e- rT E(S T - K) + . 

The parameters for the Heston model are 

v = 0.05, k = 5.07, A = 0.0457, 9 = 0.48, T = 1 

So = 100, p = r = 0.0319, p = -0.7, K = 105. 

(Since the riskfree measure is used for the valuation we have p = r.) In view of the 
above convergence result for the log-Heston model, we use the number of levels L = 
[log 2 (Te~ 1 )] and the number of repetitions iVj = f Le~ 2 T2~ e ] , £ = 0, . . . , L, for a given 
input accuracy s > 0, see [19] . 

The table below shows the empirical root mean square error 



rmsq = 




for the Multi-level estimator versus the required number of total Euler steps. The latter 
is proportional to the overall computational cost of the estimator, i.e. the number 
of used random numbers, number of function evaluations and number of arithmetic 
operations. The P^] are iid copies of the Multi-level estimator P m \ and we use M = 5-10 4 . 
The reference value p = 7.46253 was obtained by a numerical evaluation of its Fourier 
transform representation, see e.g. [3]. 

For comparison, we also provide the corresponding numerical data for the standard 
Monte Carlo estimator P st , see fill I) , for which we use the relation A 2 = T/N to match 
stepsize A and numbers of repetitions N. For the same parameters as above the empirical 
root mean square error of P s t is again estimated using M — 5 • 10 4 repetitions. 



e 


Euler steps of P m i 


rmsq cmp of P m i 


Euler steps of P st 


rmSC lcmp of p st 


2- 3 


1056 


1.369616 


512 


1.444497 


2 -4 


7168 


0.685299 


4096 


0.714207 


2" b 


43520 


0.352762 


32768 


0.357962 


2 -6 


245760 


0.181384 


262144 


0.179231 


2- 7 


1318912 


0.093485 


2097152 


0.089618 


2-8 


6815744 


0.047139 


16777216 


0.044821 



The numerical data are in good accordance with the predicted convergence behavior, 
that is 
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• for the Multi-level estimator a root mean square error of order e for a computa- 
tional cost of order e~ 2 \ log(e)| 2 

• and for the standard estimator a root mean square error of order e for a compu- 
tational cost of order e~ 3 . 

In particular halving the input accuracy leads for both estimators (approximately) to a 
halving of the empirical root mean square error. Moreover, these results illustrate nicely 
the superiority of the Multi-level estimator for small input accuracies. 



6. Summary and Outlook 

In this article we gave a survey on recent results on the convergence of numerical 
methods for stochastic differential equations in mathematical finance. The presented 
results include: 

• the pathwise convergence of general Ito- Taylor schemes for strictly positive SDEs 
with smooth but not globally Lipschitz coefficients (Section 2); 

• the construction of structure, i.e. positivity, preserving approximation schemes 
(Sections 2 and 5); 

• the strong convergence of Euler-type methods for the CEV model and the CIR 
process (Section 3); 

• the explosion of the moments of the Euler scheme for SDEs for the 3/2-model 
(Section 3); 

• the strong convergence of the drift-implicit Euler scheme for the Ait-Sahalia 
model (Section 4); 

• strong convergence rates for the approximation of the CIR and the log-Heston 
model using a drift-implicit Euler-type method (Section 5). 

However many unsettled questions are remaining: the exact strong convergence rate 
of the Euler scheme for the CEV and CIR processes, the existence or non-existence 
of moment explosions for approximation schemes of the Heston model, how to prevent 
moment explosions (if they happen) by simple modifications of the scheme etc. And 
even if these questions are answered, the question remains whether there is a 'general 
theory' for numerical methods for SDEs from mathematical finance or do these SDEs 
have to analysed one by one. So, the numerical analysis of SDEs arising in finance will 
be still an active and challenging field of research in the future. 
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